Probabilidad de quema

103039 - Gestió de Riscos en la Planificació Forestal - UdL

Author

Pere Gelabert

Published

February 12, 2025

1 Antecedentes

El incendio forestal del Pont de Vilomara, ocurrido en julio de 2022, fue uno de los más devastadores de la temporada en Catalunya, arrasando más de 1.700 hectáreas de terreno en la comarca del Bages. Este incendio se propagó con gran rapidez debido a las altas temperaturas, la sequía extrema y la sequedad del combustible, poniendo en peligro numerosas viviendas y obligando a la evacuación de cientos de personas.

Una de las zonas más afectadas fue la urbanización River Park, donde el fuego avanzó sin control y consumió decenas de casas, dejando a muchas familias sin hogar. La intensa labor de los bomberos permitió contener las llamas antes de que la destrucción fuera aún mayor, pero los daños materiales y el impacto fue significativo. River Park se convirtió en un símbolo de la vulnerabilidad de las urbanizaciones ubicadas en la Interfase Urbano Forestal ante el riesgo creciente de incendios forestales, generando un debate sobre la necesidad de mejorar las medidas de prevención y gestión del territorio.

Por ello, este trabajo se centra en la simulación de la longitud de llama, por ende la energía liberada, y en la propuesta de tratamientos que permitan reducirla en las inmediaciones de las zonas habitadas. A través de este ejercicio, buscamos identificar medidas efectivas de gestión del combustible y planificación del territorio que minimicen la propagación del fuego y protejan a las comunidades ante futuros incendios.

2 Archivo Landscape

El proyecto PREVINCAT es un proyecto liderado por el Centro de Ciencia y Tecnología Forestal de Cataluña (CTFC), junto con los Bomberos de la Generalitat, para proporcionar información básica para la modelización de incendios forestales y planificar futuras actuaciones para la prevención de incendios forestales.

Un archivo de paisaje (.LCP) consiste en un ráster en formato multi-banda, que se utiliza para las simulaciones del comportamiento y efectos de los incendios forestales. Concretamente la información que debe contener es:

  • Elevación

  • Pendiente

  • Orientación

  • Modelos de combustible (Scott and Burgan, 2005 - Tabla en Campus Virtual)

  • Cobertura arbórea

  • Altura del arbolado

  • Altura de la primera rama

  • Densidad aparente de las copas

    Variables que componen el fichero Landscape (.lcp) Fuente: Finney, 2005

2.1 Definición del área de interes (AOI)

El proyecto PREVINCAT nos permite acceder a dicha información para toda Catalunya. PREVINCAT utiliza las Zonas Homogéneas de Régimen de Incendio (ZHR) para distribuir la información. Estas ZHR son descargables en el siguiente enlace y debéis descargar las que intersecten con nuestra área de estudio: El área quemada en el IF Pont de Vilomara + un buffer de 2km.

Para ello, primero de todo cargaremos las librerías necesarias para trabajar. Recordad que si es la primera vez que utilizáis la librería debéis instalarla primero utilizando la función install.packages().

library(sf)      # Manejo y análisis de datos espaciales en formato vectorial
library(terra)   # Procesamiento y análisis de datos raster
library(gstat)   # Modelado geoestadístico e interpolación espacial
library(tidyverse)  # Conjunto de paquetes para manipulación, visualización y análisis de datos
library(tmap) # Producción cartogràfica en formato visor y plot
library(spatialEco) # Analisis espacial avanzado

2.2 Acceso y mosaicado información previncat

En segundo lugar cargaremos nuestra área de estudio, el Incendio Forestal de el Pont de Vilomara, disponible en Campus Virutal, le aplicaremos un buffer de 2km y lo definiremos como AOI. Se aplica dicha distancia para evitar el efecto borde que puede producirse en los limites del AOI.

AOI <- st_read("./Imputs/Perimetre_PontVilomara_25831.shp",quiet=T) %>% st_buffer(2000)

En tercera instancia, se debe cargar y analizar que ZHR intersectan con nuestro AOI. Para ello, primeramente se carga el shapefile de con las delimitaciones espaciales de la ZHR utilizando la función st_read() y posteriormente se evalúa la intersección con la función st_filter().

ZHR <- st_read("./Preparacion/ZHR2014/ZHR.shp", quiet=T)

ZHR_sel <- ZHR %>% st_filter(AOI)
message(paste("ZHR seleccionadas:",paste(ZHR_sel$ZHR, collapse = ",")))
ZHR seleccionadas: 26,36,37
tm_shape(ZHR)+
  tm_polygons()+
  tm_shape(ZHR_sel)+
  tm_polygons(col="red")+
  tm_shape(AOI %>% st_cast("LINESTRING")) +
  tm_lines(col="blue")

Como se puede observar el área de estudio intersecciona con las ZHR 26,36,37 y debemos decargar todos los rasters de las 6 vairables mencionadas anteriormente. Una vez descargados se procederá al mosaicado y recorte por el AOI.

Advertencia

Los archivos deben guardarse y descomprimirse en una única carpeta

#Ruta donde se almacenan los archicos descargados y descomprimidos
path <- "./Previncat/"

# Creación de vector con los nombres de las variables
vars <- list.files(path,
           pattern = ".asc$",
           full.names = F) %>% substr(.,1,str_locate(.,"_")[,1]-1) %>% 
  unique(.)

# Creación de vector con las rutas a los archivos
files <- list.files(path,
                    pattern = ".asc$",
                    full.names = T) 

#Mosaicar y recortar
for(i in vars){
  r.files <- files[files %>% str_detect(i)]
  r <- lapply(r.files,rast)
  m <- mosaic(sprc(r)) %>% crop(.,AOI)
  plot(m, main=i)
  names(m) <- i
  m <- project(m, "epsg:25831","near")
  writeRaster(m,paste0("./LCP_vars/",i,".tif"),overwrite=T)
}

3 Puntos de ignición simulados

En el presente ejercicio simularemos 10.000 incendios iniciados en los puntos con mayor probabilidad de inicarse un incendio según el histórico de incendios.

Para ello debemos:

  1. Seleccionar los puntos de ignición ocurridos en un radio de 10km del incendio

  2. Estimar el umbral de distancia a partir del cual los incendios son diferentes en términos de área quemada

  3. Calcular la densidad kernel normalizada de igniciones con una distancia de búsqueda establecida en el punto 2. Además se debe establecer como densidad 0 aquellos espacios no combustibles.

  4. Generar una muestra aleatoria balanceada espacialmente según la densidad observada en el punto 3.

Info

Encontraréis los puntos de ignición registrados en la EGIF (1998-2015) en el apartado Recursos de Campus virtual

3.1 Puntos de ignición ocurridos en un radio de 10km

Para ampliar el rango espacial de análisis y tener más muestra que las posibles igniciones ocurridas en el AOI, se amplía la búsqueda de igniciones ocurridas dentro de un radio de 10km.

buff <- AOI %>% st_buffer(8000) # Previamente habiamos aplicado 2km al general el AOI

ignit <- st_read("./Imputs/Ignitions/Ignicions.shp", quiet=T)
ignit_20 <- ignit %>% st_filter(buff)

tm_shape(ignit)+
  tm_dots("grey")+
  tm_shape(ignit_20)+
  tm_dots("brown")

#Eliminar de la memoria todas las igniciones
rm(ignit)
# Vaciado de la memoria de papelera
gc()
          used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells 3579957 191.2    6996292 373.7  6996292 373.7
Vcells 5680972  43.4   12345926  94.2 12345925  94.2

3.2 Dependencia espacial del área quemada

El área quemada es un indicador de la forma en la que queman los incendios y de las características topográficas, climatológicas y de vegetación que rigen en una zona. Es por ello que se busca ese umbral de distancia a partir del cual los incendios muestran un comportamiento diferente. Ese valor será usado como rádio de búsqueda en la función de densidad desarrollada en el siguiente punto.

Para ello, primero se calibra el semivariograma. Donde se representa la semivarianza a lo largo del gradiente de distancia. Esta, por lo general, tiende a crecer hasta que se estabiliza. El punto de estabilización nos marca ese umbral de distancia a partir del cual los incendios son diferentes en términos de área quemada y todo las variables que subyacen bajo ella. Se calcula el semivariograma empírico que es derivado de muestras observaciones, de la siguiente manera:

  • Se establece una relación o formula variable~1.

  • La fuente de información

  • cutoff hace referencia a la distancia máxima a evaluar

  • width a los cortes cada \(x\) unicades de distancia

ve <- variogram(SizeHa~1,ignit_20,cutoff=10000, width=250)
head(ve)
   np      dist       gamma dir.hor dir.ver   id
1 107  133.9027   19.560791       0       0 var1
2 142  375.0521    7.853788       0       0 var1
3 200  625.5520  395.953401       0       0 var1
4 217  877.3185  697.608888       0       0 var1
5 314 1129.0535 1508.768882       0       0 var1
6 342 1375.3032 2207.342806       0       0 var1
plot(ve, plot.numbers = T)

Una vez tenemos el semivariograma empírico, ajustamos el semivariograma teórico. Que representa una curva con una ecuación conocida. Debemos seleccionar aquella curva que se ajuste mejor a la distribución de puntos observada en el semivariograma empírico:

En este caso en particular, el semivariograma spherical (Sph) es el que ajusta mejor

vt <- fit.variogram(ve, vgm("Sph"))
plot(ve, pl = T, model = vt)

kableExtra::kable(vt%>% select(model,range))
model range
Nug 0.000
Sph 7494.745

Como se observa en la tabla anterior, la distancia a la que los incendios empiezan a presentar diferencias en aria quemada es a los ~7,5km.

3.3 Densidad de puntos

Finalmente se deben generar una densidad de puntos donde se debe establecer un rádio de búsqueda máximo de 7,5km.

# Establecer como referencia espacial un raster cualquiera del archivo LCP
r.ref <- rast("./Preparacion/LCP_vars/cbd.tif")

# Cálculo de la densidad kernel
dens <- sf.kde(
  x = ignit_20,
  ref = r.ref,
  res = 20,
  bw = vt$range[2],
  standardize = T,
  scale.factor = 1
) %>% resample(r.ref,"near")

UBmask <- rast("./Preparacion/LCP_vars/model.tif")
UBmask[UBmask<=100] <- 0
UBmask[UBmask!=0] <- 1

plot(UBmask)

dens <- dens*UBmask

names(dens) <- "dens"
varnames(dens) <- "kernel density"
crs(dens) <- crs(r.ref)
plot(dens)

plot(AOI %>% st_geometry(),add=T, border="red")

3.4 Generación de muestra aleatoria espacialmente balanceada

Creación de una muestra aleatoria balanceada por la densidad. Las especificaciones de FlamMap requieren de los puntos de ignición simulados tengan las siguientes características:

  • nombre de las columnas id, X_coord, Y_coord

  • separador ,

  • decimal .

sample <- dens %>% as.data.frame(xy=T) %>% 
  slice_sample(n=10000, weight_by = dens) %>% 
  select(-dens) %>% 
  mutate(ID=1:nrow(.) %>% as.integer()) %>% 
  rename(X_coord=x,
         Y_coord=y) %>% 
  relocate(ID,X_coord,Y_coord)
write.table(sample,"./Preparacion/pt10000.txt",sep = ",", dec=".",row.names = F)

Adjuntando el paquete: 'tidyterra'
The following object is masked from 'package:stats':

    filter

4 Simulación del riesgo de incendio en FlamMap

4.1 Generación del archivo landscape (*.lcp)

Ejecutar en FlamMap

Landscape > Create Landscape file…

Importante

Revisar unidades en las que PREVINCAT sirve los datos y establecerlas correctamente en el desplegable

4.2 Definición de características ambientales

Los vientos y la humedad del combustible son los principales factores que propician la propagación de incendios forestales. Dado que queremos evaluar una situación extrema, fijaremos el percentil 97 de las condiciones de viento y humedad observadas como umbral. Este valor representa un escenario de alto riesgo, pero no es un caso completamente atípico. En este contexto, se elige porque permite simular situaciones de peligro elevado pero realista, proporcionando una base sólida para la planificación y gestión de incendios forestales.

Rosa de los vientos Pont de Vilomara. Fuente: Meteocat

Según los datos aportados por la Agencia Meteorológica Meteocat, la velocidad de viento dominante es N

Por otro lado, debemos extrapolar el valor del P97 según la tabla de % de datos de vv que aparece en el margen inferior de la imagen de la rosa de los vientos mediante el siguiente procedimiento

tableVV <- data.frame(inf=c(0,.5,2.5,5,7.5,10,20),
           sup=c(0.5,2.5,5,7.5,10,20,Inf),
           freq_datos=c(18.8,62.4,16.3,2.5,0,0,0)) %>% 
  mutate(frec_accum=cumsum(freq_datos))

print(tableVV)
   inf  sup freq_datos frec_accum
1  0.0  0.5       18.8       18.8
2  0.5  2.5       62.4       81.2
3  2.5  5.0       16.3       97.5
4  5.0  7.5        2.5      100.0
5  7.5 10.0        0.0      100.0
6 10.0 20.0        0.0      100.0
7 20.0  Inf        0.0      100.0

\(P_k = L_i + \left( \frac{k - F_{i-1}}{F_i - F_{i-1}} \right) \times (L_s - L_i)\)

  • \(P_k\)​ valor de la variable en el percentil deseado.

  • \(L_i\)​ es el límite inferior del intervalo donde se encuentra \(P_k\)​.

  • \(L_s\)​ es el límite superior del intervalo.

  • \(F_{i-1}\) es la frecuencia acumulada antes del intervalo.

  • \(F_i\)​ es la frecuencia acumulada al final del intervalo.

  • \(k\) es el percentil deseado (en este caso, 97).

calcular_percentil <- function(tabla, percentil) {
  # Encontrar el intervalo donde se encuentra el percentil
  for (i in 1:nrow(tabla)) {
    if (tabla$frec_accum[i] >= percentil) {
      # Extraer valores del intervalo correspondiente
      Li <- tabla$inf[i]      # Límite inferior
      Ls <- tabla$sup[i]      # Límite superior
      Fi_1 <- ifelse(i == 1, 0, tabla$frec_accum[i - 1])  # Frecuencia acumulada anterior
      Fi <- tabla$frec_accum[i]  # Frecuencia acumulada actual
      
      # Aplicar la fórmula de interpolación lineal
      Pk <- Li + ((percentil - Fi_1) / (Fi - Fi_1)) * (Ls - Li)
      
      return(Pk)
    }
  }
  return(NA)  # Devuelve NA si el percentil no se encuentra en la tabla
}

VVP97 <- calcular_percentil(tableVV,97)
print(VVP97)
[1] 4.923313

Por tanto la velocidad extrema de viento se sitúa a los \(4.92 m/s \space (17.712 km/h)\) y la dirección de viento dominante es Norte.

4.2.1 Determinar niveles de humedad del combustible

Por otro lado, debe establecerse el contenido de humedad en la vegetación. Este se calcula mediante indices y ecuaciones que contemplan trayectorias temporales de humedad relativa, temperatura y precipitación.

En este caso utilizaremos los percentiles de humedad del combustible estimados para la Vegueria de la Catalunya Central por Alcasena, 2019.

Percentil 97 de humedades de combustible y de velocidad de viento
FuelModel x1h x10h x100h LW LH VV DV
0 8 9 13 50 60 4.92 N

5 Determinar la duración de los incendios simulados

Para intentar recrear al máximo la realidad de la zona de estudio, se buscará cuan es el tiempo optimo en alcanzar el tamaño medio de los incendios de la zona.

Por ello se buscan los incendios que hemos determinado anteriormente que eran parecidos en area quemada, aquellos que se encuentran a ~7,5 Km del AOI y calculamos la media de sus superficies quemadas.

ignit_20 %>% st_filter(AOI %>% st_buffer(vt$range[2])) %>% 
  st_drop_geometry() %>% 
  filter(SizeHa>10) %>% 
  summarise(mean(SizeHa)) %>% pull()
[1] 101.3135

Por tanto debemos recrear 10.000 incendios de unas que de promedio quemen 101 ha con las condiciones ambientales extremas establecidas establecidas. Para ello iremos incrementando el tiempo de quema y hasta alcanzar el umbral de tiempo necesario para quemar esa misma área promedio.

Ejemplo de otra area de estudio donde se pretendia alcanzar una área quemada de 293 ha y se alcanzó con 303min

Ejecución algoritmo MTT
Ejecutar en FlamMap

Analysis Area > New FlamMap/MTT/TOM run con los parámetros

Debéis ajustar todos los parámetros requeridos (marcados en rojo)

Finalmente, debéis repitiendo simulaciones e ir incrementando el tiempo hasta alcanzar una área media quemada similar al histórico.

Hay ejecutaremos el algoritmo Minimum Travel Time (MTT), que, de manera simplificada, calcula el comportamiento del fuego buscando el conjunto de caminos con tiempos mínimos de propagación del fuego a partir de igniciones proporcionadas (ampliar información). Estableciendo el tiempo de quema calibrado y marcando las casillas para seleccionar los outputs deseados

Outputs a seleccionar
  • Burn probabilities

  • Perímeters

  • FLP-20 Bin Metric

  • Fire Size List

Muy importante

Se deben exportar los archivos des de el panel lateral de FlamMap y abrirlos con QGIS o ArcGIS, para los objetos espaciales y Excel o R para las tablas. No se recomienda abrirlos en el propio software ya que al ser pesado el programa puede no responder y perderse todo el trabajo realizado.


Parte 2: En esta segunda parte se procesan, tratan e interpretan los resultados obtenidos.

6 Outputs MTT

En el siguiente punto, se listan los productos Los productos que nos devuelve el algoritmo MTT:

  • Perimetros MTT: Los 10000 perímetros de las diferentes igniciones simuladas

  • Burn Probabilities: Ráster con la probabilidad

  • FLP metric: Longitud de llama condicional. Probabilidad de tener una longitud de llama determinada.

  • Fire Size List: Tabla con las coordenadas de cada ignición y el tamaño de incendio en cada una de ellas.
FIRE_NUM XStart YStart FireSize.Acres FireSize.Ha
0 410595.5 4621331 362.9484 146.88013
1 409436.7 4613521 236.7270 95.80008
2 410649.4 4621957 537.8002 217.64019
3 410258.6 4624639 594.9309 240.76021
4 406216.4 4614296 196.2017 79.40007
5 404922.9 4616681 323.5104 130.92012

La longitud de llama responde a parámetros físicos como la intensidad de la llama (kW/m) que está ligada al calor de combustión o a la tasa de consumo de la vegetación. Se estima a partir de la expresión propuesta por Byram (1959):

\[ L_f=0.0775·I_b^{0.46} \]

Donde:

  • \(L_f\) es la longitud de llama

  • \(I_b\) es la intensidad de Byram

6.1 Cálculo longitud de llama

# Leer los datos de Probabilidad de Longitud de Llama (FLP) y ordenarlos en orden descendente por PBurn
LFP <- read_csv("./Outputs/FLP.csv") %>% arrange(desc(PBurn))
Rows: 256055 Columns: 23
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (23): XPos, YPos, PBurn, FIL1, FIL2, FIL3, FIL4, FIL5, FIL6, FIL7, FIL8,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Crear una secuencia de constantes para el cálculo de la longitud de llama
constantes <- c(seq(0.5,9.5,.5)-0.25, 12)

# Calcular la longitud de llama (FL) y convertir los datos a un objeto espacial 'sf'
FL <- LFP %>% select(XPos,YPos,PBurn) %>% bind_cols(
  sweep(LFP[, 4:ncol(LFP)], 2, constantes, `*`) %>% mutate(FL=rowSums(.)) %>% select(FL)) %>% 
  st_as_sf(coords = c("XPos","YPos")) %>% 
  st_set_crs(25831)  # Establecer el sistema de coordenadas

# Leer el raster de referencia y ajustar su sistema de coordenadas
ras.ref <- r.ref
crs(ras.ref) <- "epsg:25831"

# Rasterizar los datos de longitud de llama (FL) utilizando el raster de referencia
FL <- rasterize(FL, ras.ref, field="FL") %>% mask(.,AOI)

# Visualizar el raster resultante
plot(FL %>% crop(AOI))

class_mat <- matrix(c(-Inf,1,1,
                      1,2.5,2,
                      2.5,3.5,3,
                      3.5,Inf,4),
                    ncol=3, byrow = T)
FL_clas <- classify(FL,class_mat)



coltab(FL_clas) <- data.frame(values = c(1,2,3,4),
                        cols = c("darkgreen", "gold", "orange", "brown"))

mylevels <- 

levels(FL_clas) <- data.frame(values = c(1,2,3,4),
                       cats = c("Bajo", "Medio", "Alto", "Muy Alto"))

plot(FL_clas)

RiverPark<- st_read("./Imputs/RiverPark.shp", quiet=T)
RiverPark_fill <- st_multipolygon(lapply(st_geometry(RiverPark) , function(x) x[1]))
plot(FL_clas %>% crop(RiverPark %>% st_buffer(200)))


plot(RiverPark_fill %>% st_geometry(), , add=T, border="blue", lwd=2)

6.2 Tratamientos

Para reducir la exposición de bienes a altas intensidades de fuego se suelen realizar una serie de tratamientos forestales, como por ejemplo:

Reducción de combustibles de superficie

  • Quemas prescritas. Bosques de coníferas y pastos arbustivos. Requiere tecnificación,  riesgo de escape. Ventana meteorológica. Plan de quemas.

  • Tratamientos mecánicos. Interfaz urbano-forestal, bosques de frondosas y zonas de elevada carga de combustibles. No requiere tecnificación.

Tratamiento de fracción arbolada

  • Aplicación de poda baja.

  • Eliminación de la tangencia en copas.

Conservación de áreas tratadas

  • Uso de ganadería extensiva en sierras o montes cercados. Grandes superficies.

  • Uso de pastoreo dirigido para la conservación de cortafuegos o puntos estratégicos

En este caso concreto estos tratamientos se podrían plasmar como:

  1. Alterar el Modelo de combustible, por ejemplo de SB4 a SB1

  2. Reducir la FCC hasta el 50%

  3. Subir la base de la copa hasta \(3 m\)

  4. Reducir la densidad aparente de copas hasta \(0.08 kg/m^3\)

  5. Reducir la altura del arbolado mediante claras selectivas.

7 Ejercicio a realizar

  1. Evaluar cual es la longitud de llama en las proximidades de la urbanización (un buffer de 60 m 3 píxeles).

  2. Explorar cuales son los valores de las variables utilizadas en el archivo *.lcp en las inmediaciones de la urbanización (buffer 60m).

  3. Diseñar tratamientos para reducir la longitud de llama en las inmediaciones de la urbanización.

  4. Replantear los rasters del archivo landscape y reflejando los tratamientos.

  5. Correr de nuevo la simulación de 10000 igniciones con FamMap incorporando los nuevos rásters